Spline capa delgada

library(fields) # esta librería es la que nos ayudará a ajustar el spline.
library(rgl)    # esta librería nos permite graficar en 3D.

# Datos, ahí nada más hay que moverle a la ruta
data <- read.csv("/Users/fernandoperezmillan/Desktop/Aprendizaje Estadístico/Punto extra/Figura+1+8a+Datos.csv")

# Extraemos las columnas
x <- data$Education
y <- data$Seniority
z <- data$Income

# Generamos las x's y y's sobre las cuales vamos a interpolar
grid_x <- seq(min(x), max(x), length = 100)
grid_y <- seq(min(y), max(y), length = 100)

# Interpolación usando un Rough Thin Plate Spline (lo hacemos "rough" escogiendo un lambda cercano a 0 o 0).
# Usamos la librería fields, en particular la función Tps.
lambda_value <- 0
interpolation <- Tps(cbind(x, y), z, lambda = lambda_value)
z_interp <- predictSurface(interpolation, list(grid_x, grid_y))

# Abrimos el graficador 3D (para los que tienen mac, hay que instalar XQUARTZ)
#open3d()

# Graficamos los datos originales
plot3d(x, y, z, col = 'red', size = 0.7, type = "s",
       xlab = "Years of Education", ylab = "Seniority", zlab = "Income",
       xlim = range(grid_x), ylim = range(grid_y), zlim = range(z))


# Añadimos la superficie que interpolamos
surface3d(z_interp$x, z_interp$y, z_interp$z, col = terrain.colors(100), alpha = 0.9)

# Añadimos un gris azul sobre la superficie que interpolamos, para que se vea más mono
for (i in 1:length(grid_x)) {
  lines3d(rep(grid_x[i], length(grid_y)), grid_y, z_interp$z[i, ], col = rgb(0, 0, 1, 0.3), lwd = 1)
}
for (j in 1:length(grid_y)) {
  lines3d(grid_x, rep(grid_y[j], length(grid_x)), z_interp$z[, j], col = rgb(0, 0, 1, 0.3), lwd = 1)
}

# Ajustamos el ángulo de la gráfica para que se vea mejor
# view3d(theta = 45, phi = 30, zoom = 0.8)
view3d(userMatrix = rotationMatrix(pi/4, 1, 0, 0) %*% rotationMatrix(pi/3, 0, 1, 0))

# Añadimos un grid en los ejes
grid3d(c("x", "y", "z"))

# Otros ajustes
axes3d(edges = "bbox")
bg3d("white")
rgl.bringtotop()

# Embed the RGL plot into the HTML using rglwidget()
rglwidget()

#### LECTURA DE DATOS Y LIBRERÍAS ####
# ********************

# Aquí cargamos las librerías necesarias para nuestro análisis
library(fields)
library(ggplot2)
library(gridExtra)
library(grid)
library(spam)
library(dplyr)

# Leemos los datos desde el archivo CSV
datos <- read.csv("/Users/fernandoperezmillan/Desktop/Figura 1 8a Datos.csv", header = TRUE, fileEncoding="latin1")
head(datos)

#### PANEL (a): Thin-Plate Spline (TPS) ####
# *************************************************************

# La figura 2.6 del ISLR indica que ajusta un rough spline (zero fit error on the training data).
# Por eso, aquí ajustamos el modelo TPS sin penalización (lambda = 0).
# Se tomó la función Tps porque fue la que mejor funcionó para ajustar los datos.
# Nos va a permitir crear la superficie suave que queremos.

modelo_tps <- Tps(cbind(datos$Education, datos$Seniority), datos$Income, lambda = 0)

# Generamos la cuadrícula para la predicción utilizando los rangos de educación y antigüedad

rango_educacion <- range(datos$Education)
grid_educacion <- seq(rango_educacion[1], rango_educacion[2], length.out = 30)
rango_antiguedad <- range(datos$Seniority)
grid_antiguedad <- seq(rango_antiguedad[1], rango_antiguedad[2], length.out = 30)
cuadricula <- expand.grid(Educacion = grid_educacion, Antiguedad = grid_antiguedad)

# Predecimos los valores de ingreso en la cuadrícula generada
prediccion_ingreso <- matrix(predict(modelo_tps, cuadricula), 30, 30)

# Graficamos el modelo TPS con las predicciones sobre la cuadrícula
grafico_tps <- persp(grid_educacion, grid_antiguedad, prediccion_ingreso, theta = 30, phi = 30,
                     col = "gold",  # Color amarillo brillante para replicar la figura 2.6 del ISLR
                     border = "darkblue",  # Líneas de cuadrícula azules para replicar la figura 2.6 del ISLR
                     lwd = 0.2,  # Líneas de cuadrícula delgadas para replicar la figura 2.6 del ISLR
                     expand = 0.6, shade = 0.4,  # Ajuste para una apariencia más suave como la de la figura 2.6 del ISLR
                     xlab="Años de estudio", ylab="Experiencia", zlab="Ingreso anual",
                     main="Ajuste de un modelo no paramétrico: \n Spline de Capa Delgada (SCD)")

# Aquí graficamos los puntos originales sobre la superficie
# Agregamos los puntos a la gráfica previamente hecha y le damos formato a los puntos
puntos_obs <- trans3d(datos$Education, datos$Seniority, datos$Income, grafico_tps) 
points(puntos_obs, col = "red", pch = 16)

# Ahora graficamos las predicciones del modelo conectando con segmentos las predicciones y los valores reales
puntos_pred <- trans3d(datos$Education, datos$Seniority, fitted(modelo_tps), grafico_tps)
segments(puntos_obs$x, puntos_obs$y, puntos_pred$x, puntos_pred$y)